Note: Open scripts.Rproj first, then script. To easily use relative paths, click the down button next to knit and then click “Knit Directory –> Project Directory”. This should make loading and saving files much easier.

WGCNA Analysis

This script is based off of Langfelder P, Horvath S (2008) WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 2008, 9:559 (link to paper) Red text indicates regions that require the user to modify. Regardless, the user should check over all code blocks to ensure that everything is running correctly.

1. Load packages

Load packages

library(tidyverse, quietly = TRUE)
library(genefilter, quietly = TRUE)
library(edgeR, quietly = TRUE)
library(RColorBrewer, quietly = TRUE)
library(WGCNA, quietly = TRUE)
library(flashClust, quietly = TRUE)
library(gridExtra, quietly = TRUE)
library(ComplexHeatmap, quietly = TRUE)
library(goseq, quietly = TRUE)
library(dplyr, quietly = TRUE)
library(clusterProfiler, quietly = TRUE)
library(simplifyEnrichment, quietly = TRUE)

2. Set variables

Change file names and conditions where appropriate.

# Input unfiltered data
treatmentinfo.file <- "samples.tsv"
gcount.file <- "salmon.numreads.tsv"
# Output DEG results
WGCNA_results.file <- "salmon.numreads.WGCNA_results.tsv"
WGCNA_samples.file <- "salmon.numreads.WGCNA_samples.tsv"

# pOverA gene filtering thresholds (only process genes with over these cutoffs)
pOverA.cutoff.P <- 0.8 # >80% of samples
pOverA.cutoff.A <- 10  # >10 read counts

# Column name to use for analysis of data using WGCNA
group <- "group"

# Minimum module size - I chose 30 as it is the default value chosen by most studies using WGCNA.
minModuleSize <- 30

3. Load, clean, and pre-processing datasets

Load the input file containing the treatment information.

treatmentinfo <- read.csv(treatmentinfo.file, header = TRUE, sep = "\t", fileEncoding="UTF-8-BOM") # Read in file
head(treatmentinfo)
##              sample_id unit lib_type         fq1         fq2 species plug_id
## 1 Pacuta_ATAC_TP1_1043    1       pe SRR14610932 SRR14610932  Pacuta    1043
## 2 Pacuta_ATAC_TP1_1775    1       pe SRR14610931 SRR14610931  Pacuta    1775
## 3 Pacuta_ATAC_TP1_2363    1       pe SRR14610930 SRR14610930  Pacuta    2363
## 4 Pacuta_ATAC_TP3_1041    1       pe SRR14610929 SRR14610929  Pacuta    1041
## 5 Pacuta_ATAC_TP3_1471    1       pe SRR14610928 SRR14610928  Pacuta    1471
## 6 Pacuta_ATAC_TP3_1637    1       pe SRR14610927 SRR14610927  Pacuta    1637
##   treatment timepoint            reef ploidy  group
## 1      ATAC         1 Lilipuna.Fringe      3 Group4
## 2      ATAC         1      Reef.42.43      2 Group5
## 3      ATAC         1         Reef.18      3 Group3
## 4      ATAC         3      Reef.11.13      3 Group2
## 5      ATAC         3      Reef.35.36      3 Group3
## 6      ATAC         3      Reef.11.13      3 Group2
# Check we have the right column names
headers <- c("sample_id")
if( all(headers %in% colnames(treatmentinfo)) ){
  cat(paste(treatmentinfo.file, "has the required columns:", paste(headers, collapse=', ')))
} else {
  stop(paste(treatmentinfo.file, "is missing required columns:", paste(headers, collapse=', ')))
}
## samples.tsv has the required columns: sample_id

Filter treatmentinfo to just the samples that we want to analyze. Generally we would use all samples for WGCNA, however, in some cases we might want to split the data if we have samples from different experiments. Adjust filtering as required.

#treatmentinfo <- treatmentinfo %>%
#  filter(condition %in% c("ATAC", "HTAC") & timepoint=="8") %>%
#  mutate(condition = factor(condition, levels = c("ATAC", "HTAC")))
head(treatmentinfo)
##              sample_id unit lib_type         fq1         fq2 species plug_id
## 1 Pacuta_ATAC_TP1_1043    1       pe SRR14610932 SRR14610932  Pacuta    1043
## 2 Pacuta_ATAC_TP1_1775    1       pe SRR14610931 SRR14610931  Pacuta    1775
## 3 Pacuta_ATAC_TP1_2363    1       pe SRR14610930 SRR14610930  Pacuta    2363
## 4 Pacuta_ATAC_TP3_1041    1       pe SRR14610929 SRR14610929  Pacuta    1041
## 5 Pacuta_ATAC_TP3_1471    1       pe SRR14610928 SRR14610928  Pacuta    1471
## 6 Pacuta_ATAC_TP3_1637    1       pe SRR14610927 SRR14610927  Pacuta    1637
##   treatment timepoint            reef ploidy  group
## 1      ATAC         1 Lilipuna.Fringe      3 Group4
## 2      ATAC         1      Reef.42.43      2 Group5
## 3      ATAC         1         Reef.18      3 Group3
## 4      ATAC         3      Reef.11.13      3 Group2
## 5      ATAC         3      Reef.35.36      3 Group3
## 6      ATAC         3      Reef.11.13      3 Group2

Load the input file containing the gene count matrix.

gcount <- as.data.frame(read.csv(gcount.file, header = TRUE, sep = "\t", fileEncoding="UTF-8-BOM")) # Read in file

# Check we have the right column names
headers <- c("Name", treatmentinfo$sample_id)
if( all(headers %in% colnames(gcount)) ){
  cat(paste(gcount.file, "has the required columns:", paste(headers, collapse=', ')))
} else {
  stop(paste(gcount.file, "is missing required columns:", paste(headers, collapse=', ')))
}
## salmon.numreads.tsv has the required columns: Name, Pacuta_ATAC_TP1_1043, Pacuta_ATAC_TP1_1775, Pacuta_ATAC_TP1_2363, Pacuta_ATAC_TP3_1041, Pacuta_ATAC_TP3_1471, Pacuta_ATAC_TP3_1637, Pacuta_ATAC_TP4_1060, Pacuta_ATAC_TP4_1762, Pacuta_ATAC_TP4_2002, Pacuta_ATAC_TP5_1059, Pacuta_ATAC_TP5_1563, Pacuta_ATAC_TP5_1757, Pacuta_ATAC_TP6_1050, Pacuta_ATAC_TP6_1468, Pacuta_ATAC_TP6_1542, Pacuta_ATAC_TP7_1047, Pacuta_ATAC_TP7_1445, Pacuta_ATAC_TP7_2413, Pacuta_ATAC_TP8_1051, Pacuta_ATAC_TP8_1755, Pacuta_ATAC_TP8_2012, Pacuta_ATAC_TP9_1141, Pacuta_ATAC_TP9_1594, Pacuta_ATAC_TP9_2357, Pacuta_ATAC_TP10_1159, Pacuta_ATAC_TP10_1559, Pacuta_ATAC_TP10_1641, Pacuta_ATAC_TP11_1103, Pacuta_ATAC_TP11_1777, Pacuta_ATAC_TP11_2306, Pacuta_ATHC_TP1_1207, Pacuta_ATHC_TP1_2743, Pacuta_ATHC_TP1_2977, Pacuta_ATHC_TP3_1219, Pacuta_ATHC_TP3_2534, Pacuta_ATHC_TP3_2750, Pacuta_ATHC_TP4_1220, Pacuta_ATHC_TP4_2733, Pacuta_ATHC_TP4_2993, Pacuta_ATHC_TP5_1296, Pacuta_ATHC_TP5_2212, Pacuta_ATHC_TP5_2877, Pacuta_ATHC_TP6_1254, Pacuta_ATHC_TP6_2870, Pacuta_ATHC_TP6_2999, Pacuta_ATHC_TP7_1281, Pacuta_ATHC_TP7_2409, Pacuta_ATHC_TP7_2878, Pacuta_ATHC_TP8_1459, Pacuta_ATHC_TP8_2564, Pacuta_ATHC_TP8_2861, Pacuta_ATHC_TP9_1451, Pacuta_ATHC_TP9_2873, Pacuta_ATHC_TP9_2979, Pacuta_ATHC_TP10_1205, Pacuta_ATHC_TP10_2197, Pacuta_ATHC_TP10_2550, Pacuta_ATHC_TP11_1147, Pacuta_ATHC_TP11_2668, Pacuta_ATHC_TP11_2879, Pacuta_HTAC_TP1_1653, Pacuta_HTAC_TP1_2005, Pacuta_HTAC_TP1_2414, Pacuta_HTAC_TP3_1617, Pacuta_HTAC_TP3_1642, Pacuta_HTAC_TP3_2026, Pacuta_HTAC_TP4_1581, Pacuta_HTAC_TP4_1701, Pacuta_HTAC_TP4_1767, Pacuta_HTAC_TP5_1303, Pacuta_HTAC_TP5_1571, Pacuta_HTAC_TP5_1707, Pacuta_HTAC_TP6_1330, Pacuta_HTAC_TP6_1466, Pacuta_HTAC_TP6_1744, Pacuta_HTAC_TP7_1487, Pacuta_HTAC_TP7_1728, Pacuta_HTAC_TP7_2072, Pacuta_HTAC_TP8_1329, Pacuta_HTAC_TP8_1765, Pacuta_HTAC_TP8_2513, Pacuta_HTAC_TP9_1302, Pacuta_HTAC_TP9_1486, Pacuta_HTAC_TP9_1696, Pacuta_HTAC_TP10_1225, Pacuta_HTAC_TP10_1536, Pacuta_HTAC_TP10_2064, Pacuta_HTAC_TP11_1582, Pacuta_HTAC_TP11_1596, Pacuta_HTAC_TP11_1647, Pacuta_HTHC_TP1_1239, Pacuta_HTHC_TP1_1676, Pacuta_HTHC_TP1_2210, Pacuta_HTHC_TP3_1227, Pacuta_HTHC_TP3_1418, Pacuta_HTHC_TP3_2527, Pacuta_HTHC_TP4_1169, Pacuta_HTHC_TP4_1343, Pacuta_HTHC_TP4_2195, Pacuta_HTHC_TP5_1168, Pacuta_HTHC_TP5_1415, Pacuta_HTHC_TP5_2087, Pacuta_HTHC_TP6_1138, Pacuta_HTHC_TP6_1595, Pacuta_HTHC_TP6_1721, Pacuta_HTHC_TP7_1090, Pacuta_HTHC_TP7_1427, Pacuta_HTHC_TP7_1820, Pacuta_HTHC_TP8_1184, Pacuta_HTHC_TP8_1709, Pacuta_HTHC_TP8_2304, Pacuta_HTHC_TP9_1131, Pacuta_HTHC_TP9_2202, Pacuta_HTHC_TP9_2305, Pacuta_HTHC_TP10_1238, Pacuta_HTHC_TP10_1732, Pacuta_HTHC_TP10_2300, Pacuta_HTHC_TP11_1416, Pacuta_HTHC_TP11_2185
# Cleanup gcount data
gcount <- gcount %>% 
  column_to_rownames("Name") %>% # Makes "Name" column rownames
  round() %>% # Round
  select(treatmentinfo$sample_id) # Select just columns in filtered treatmentinfo file

# View dataset attributes
d <- dim(gcount); print(paste("rows:",d[[1]]," columns:",d[[2]], sep=''))
## [1] "rows:33259 columns:119"
head(gcount)[,1:3]
##                                             Pacuta_ATAC_TP1_1043
## Pocillopora_acuta_KBHIv2___RNAseq.g16115.t1                   78
## Pocillopora_acuta_KBHIv2___TS.g16868.t1                      216
## Pocillopora_acuta_KBHIv2___RNAseq.g25384.t1                    4
## Pocillopora_acuta_KBHIv2___RNAseq.g30122.t1                    0
## Pocillopora_acuta_KBHIv2___RNAseq.g12817.t1                    1
## Pocillopora_acuta_KBHIv2___TS.g9097.t1                        90
##                                             Pacuta_ATAC_TP1_1775
## Pocillopora_acuta_KBHIv2___RNAseq.g16115.t1                   47
## Pocillopora_acuta_KBHIv2___TS.g16868.t1                      121
## Pocillopora_acuta_KBHIv2___RNAseq.g25384.t1                    0
## Pocillopora_acuta_KBHIv2___RNAseq.g30122.t1                    0
## Pocillopora_acuta_KBHIv2___RNAseq.g12817.t1                    0
## Pocillopora_acuta_KBHIv2___TS.g9097.t1                        12
##                                             Pacuta_ATAC_TP1_2363
## Pocillopora_acuta_KBHIv2___RNAseq.g16115.t1                   32
## Pocillopora_acuta_KBHIv2___TS.g16868.t1                      211
## Pocillopora_acuta_KBHIv2___RNAseq.g25384.t1                    3
## Pocillopora_acuta_KBHIv2___RNAseq.g30122.t1                    0
## Pocillopora_acuta_KBHIv2___RNAseq.g12817.t1                    2
## Pocillopora_acuta_KBHIv2___TS.g9097.t1                       107

Filter gcount to just genes with enough data across columns (exclude genes which are mostly zeros)

# Create filter for the counts data
filt <- filterfun(pOverA(pOverA.cutoff.P, pOverA.cutoff.A))
gfilt <- genefilter(gcount, filt)
# Identify genes to keep by count filter
keep <- gcount[gfilt,]
# Identify gene lists
keep <- rownames(keep)
# Gene count data filtered in PoverA, P percent of the samples have counts over A
gcount_filt <- as.data.frame(gcount[which(rownames(gcount) %in% keep),])
d <- dim(gcount);      print(paste("gcount      - rows:",d[[1]]," columns:",d[[2]], sep='')) # Before filtering
## [1] "gcount      - rows:33259 columns:119"
d <- dim(gcount_filt); print(paste("gcount_filt - rows:",d[[1]]," columns:",d[[2]], sep='')) # After filtering
## [1] "gcount_filt - rows:15467 columns:119"

Determine library size.

libSize.df <- data.frame(libSize=colSums(gcount_filt))

Make DGE object.

DGEdat <- DGEList(counts=as.matrix(gcount_filt),
                  samples=treatmentinfo,
                  group=treatmentinfo[[group]])
dim(DGEdat$counts)
## [1] 15467   119

4. Data normalization

DGEdat <- calcNormFactors(DGEdat)
head(DGEdat$samples)
##                       group lib.size norm.factors            sample_id unit
## Pacuta_ATAC_TP1_1043 Group4  4128082    0.9769418 Pacuta_ATAC_TP1_1043    1
## Pacuta_ATAC_TP1_1775 Group5  4557568    0.9699660 Pacuta_ATAC_TP1_1775    1
## Pacuta_ATAC_TP1_2363 Group3  4629399    0.9728287 Pacuta_ATAC_TP1_2363    1
## Pacuta_ATAC_TP3_1041 Group2  4870724    0.9988819 Pacuta_ATAC_TP3_1041    1
## Pacuta_ATAC_TP3_1471 Group3  3699133    0.9679178 Pacuta_ATAC_TP3_1471    1
## Pacuta_ATAC_TP3_1637 Group2  3459979    0.9964851 Pacuta_ATAC_TP3_1637    1
##                      lib_type         fq1         fq2 species plug_id treatment
## Pacuta_ATAC_TP1_1043       pe SRR14610932 SRR14610932  Pacuta    1043      ATAC
## Pacuta_ATAC_TP1_1775       pe SRR14610931 SRR14610931  Pacuta    1775      ATAC
## Pacuta_ATAC_TP1_2363       pe SRR14610930 SRR14610930  Pacuta    2363      ATAC
## Pacuta_ATAC_TP3_1041       pe SRR14610929 SRR14610929  Pacuta    1041      ATAC
## Pacuta_ATAC_TP3_1471       pe SRR14610928 SRR14610928  Pacuta    1471      ATAC
## Pacuta_ATAC_TP3_1637       pe SRR14610927 SRR14610927  Pacuta    1637      ATAC
##                      timepoint            reef ploidy group.1
## Pacuta_ATAC_TP1_1043         1 Lilipuna.Fringe      3  Group4
## Pacuta_ATAC_TP1_1775         1      Reef.42.43      2  Group5
## Pacuta_ATAC_TP1_2363         1         Reef.18      3  Group3
## Pacuta_ATAC_TP3_1041         3      Reef.11.13      3  Group2
## Pacuta_ATAC_TP3_1471         3      Reef.35.36      3  Group3
## Pacuta_ATAC_TP3_1637         3      Reef.11.13      3  Group2

5. Plot global gene expression

Log transform the counts matrix for the next plots

DGEdat.cpm <- DGEdat # Make a copy the edgeR dataset
DGEdat.cpm$counts <- cpm(DGEdat.cpm$counts, log=TRUE, prior.count=5) # Log transform the copy for the next plots

Plot a heatmap of sample-to-sample distances

gsampleDists_dev <- dist(t(DGEdat.cpm$counts)) # Calculate distance matrix
gsampleDistMatrix_dev <- as.matrix(gsampleDists_dev) # Distance matrix
rownames(gsampleDistMatrix_dev) <- colnames(DGEdat.cpm) # Assign row names
colnames(gsampleDistMatrix_dev) <- NULL # Assign col names
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255) # Assign colors
pheatmap(gsampleDistMatrix_dev, # Plot matrix
         clustering_distance_rows=gsampleDists_dev, # Cluster rows
         clustering_distance_cols=gsampleDists_dev, # Cluster columns
         col=colors) # Set colors

Principal component plot of samples

pca <- prcomp(t(DGEdat.cpm$counts)) # Calculate eigengenes
pc.data <- summary(out<-prcomp(t(DGEdat.cpm$counts))); pc.data
## Importance of components:
##                            PC1     PC2      PC3      PC4      PC5      PC6
## Standard deviation     27.0984 24.7536 23.00132 20.99766 16.05678 14.56444
## Proportion of Variance  0.1363  0.1137  0.09818  0.08182  0.04785  0.03937
## Cumulative Proportion   0.1363  0.2500  0.34817  0.43000  0.47784  0.51721
##                             PC7      PC8      PC9     PC10    PC11     PC12
## Standard deviation     13.69539 12.15218 11.81621 10.84512 10.4334 10.37330
## Proportion of Variance  0.03481  0.02741  0.02591  0.02183  0.0202  0.01997
## Cumulative Proportion   0.55202  0.57942  0.60534  0.62716  0.6474  0.66733
##                           PC13    PC14    PC15    PC16    PC17    PC18   PC19
## Standard deviation     9.56842 9.03839 8.90666 8.68261 8.11781 7.84050 7.1939
## Proportion of Variance 0.01699 0.01516 0.01472 0.01399 0.01223 0.01141 0.0096
## Cumulative Proportion  0.68433 0.69949 0.71421 0.72820 0.74043 0.75184 0.7614
##                           PC20    PC21    PC22    PC23   PC24    PC25   PC26
## Standard deviation     7.03702 6.63725 6.43582 6.07434 6.0078 5.89786 5.6402
## Proportion of Variance 0.00919 0.00818 0.00769 0.00685 0.0067 0.00646 0.0059
## Cumulative Proportion  0.77063 0.77881 0.78649 0.79334 0.8000 0.80649 0.8124
##                           PC27    PC28    PC29    PC30    PC31    PC32   PC33
## Standard deviation     5.47977 5.41242 5.32948 5.14454 5.00973 4.82876 4.6441
## Proportion of Variance 0.00557 0.00544 0.00527 0.00491 0.00466 0.00433 0.0040
## Cumulative Proportion  0.81797 0.82341 0.82868 0.83359 0.83825 0.84257 0.8466
##                           PC34    PC35    PC36   PC37    PC38    PC39    PC40
## Standard deviation     4.58795 4.49855 4.38008 4.3454 4.22217 4.13238 4.10966
## Proportion of Variance 0.00391 0.00376 0.00356 0.0035 0.00331 0.00317 0.00313
## Cumulative Proportion  0.85048 0.85424 0.85780 0.8613 0.86461 0.86778 0.87092
##                           PC41    PC42    PC43    PC44    PC45    PC46    PC47
## Standard deviation     4.08065 4.05213 3.90462 3.83247 3.79697 3.73699 3.72268
## Proportion of Variance 0.00309 0.00305 0.00283 0.00273 0.00268 0.00259 0.00257
## Cumulative Proportion  0.87401 0.87705 0.87988 0.88261 0.88528 0.88788 0.89045
##                           PC48    PC49    PC50    PC51    PC52    PC53    PC54
## Standard deviation     3.64532 3.62479 3.56088 3.54082 3.51107 3.44885 3.42440
## Proportion of Variance 0.00247 0.00244 0.00235 0.00233 0.00229 0.00221 0.00218
## Cumulative Proportion  0.89291 0.89535 0.89771 0.90003 0.90232 0.90453 0.90670
##                           PC55    PC56    PC57    PC58    PC59    PC60    PC61
## Standard deviation     3.37977 3.37074 3.33918 3.31446 3.26671 3.24470 3.22975
## Proportion of Variance 0.00212 0.00211 0.00207 0.00204 0.00198 0.00195 0.00194
## Cumulative Proportion  0.90882 0.91093 0.91300 0.91504 0.91702 0.91897 0.92091
##                           PC62   PC63    PC64    PC65    PC66   PC67    PC68
## Standard deviation     3.22092 3.2022 3.17807 3.16082 3.13362 3.1168 3.09767
## Proportion of Variance 0.00193 0.0019 0.00187 0.00185 0.00182 0.0018 0.00178
## Cumulative Proportion  0.92284 0.9247 0.92661 0.92847 0.93029 0.9321 0.93387
##                           PC69    PC70   PC71    PC72    PC73    PC74    PC75
## Standard deviation     3.07047 3.04402 3.0250 3.00161 2.98065 2.97374 2.96028
## Proportion of Variance 0.00175 0.00172 0.0017 0.00167 0.00165 0.00164 0.00163
## Cumulative Proportion  0.93562 0.93734 0.9390 0.94071 0.94236 0.94400 0.94563
##                           PC76   PC77    PC78    PC79    PC80    PC81    PC82
## Standard deviation     2.95028 2.9346 2.92463 2.89990 2.89350 2.86431 2.83641
## Proportion of Variance 0.00162 0.0016 0.00159 0.00156 0.00155 0.00152 0.00149
## Cumulative Proportion  0.94724 0.9488 0.95043 0.95199 0.95354 0.95507 0.95656
##                           PC83    PC84    PC85    PC86    PC87    PC88    PC89
## Standard deviation     2.81509 2.80595 2.79664 2.76909 2.76666 2.75655 2.73046
## Proportion of Variance 0.00147 0.00146 0.00145 0.00142 0.00142 0.00141 0.00138
## Cumulative Proportion  0.95803 0.95949 0.96094 0.96237 0.96379 0.96520 0.96658
##                           PC90    PC91    PC92    PC93    PC94   PC95    PC96
## Standard deviation     2.71741 2.70707 2.69606 2.69002 2.67369 2.6500 2.63861
## Proportion of Variance 0.00137 0.00136 0.00135 0.00134 0.00133 0.0013 0.00129
## Cumulative Proportion  0.96795 0.96931 0.97066 0.97200 0.97333 0.9746 0.97592
##                           PC97    PC98    PC99   PC100   PC101  PC102   PC103
## Standard deviation     2.63269 2.61371 2.60259 2.58884 2.57400 2.5481 2.52325
## Proportion of Variance 0.00129 0.00127 0.00126 0.00124 0.00123 0.0012 0.00118
## Cumulative Proportion  0.97721 0.97848 0.97974 0.98098 0.98221 0.9834 0.98459
##                          PC104   PC105   PC106   PC107   PC108   PC109   PC110
## Standard deviation     2.50312 2.49367 2.46185 2.44319 2.41778 2.40138 2.38468
## Proportion of Variance 0.00116 0.00115 0.00112 0.00111 0.00108 0.00107 0.00106
## Cumulative Proportion  0.98576 0.98691 0.98804 0.98914 0.99023 0.99130 0.99235
##                          PC111   PC112   PC113   PC114   PC115   PC116   PC117
## Standard deviation     2.36155 2.33573 2.32859 2.28930 2.26562 2.22176 2.21337
## Proportion of Variance 0.00103 0.00101 0.00101 0.00097 0.00095 0.00092 0.00091
## Cumulative Proportion  0.99339 0.99440 0.99541 0.99638 0.99733 0.99825 0.99916
##                          PC118    PC119
## Standard deviation     2.12894 1.91e-13
## Proportion of Variance 0.00084 0.00e+00
## Cumulative Proportion  1.00000 1.00e+00
plot(out)

#DGEdat.cpm$samples$timepoint <- gsub("TP", "", DGEdat.cpm$samples$timepoint)
DGEdat_PCcor <- lapply(DGEdat.cpm$samples, as.factor)
#DGEdat_PCcor <- as.tibble(lapply(DGEdat_PCcor, as.numeric))
DGEdat_PCcor <- as.tibble(DGEdat_PCcor)
## Warning: `as.tibble()` was deprecated in tibble 2.0.0.
## ℹ Please use `as_tibble()` instead.
## ℹ The signature and semantics have changed, see `?as_tibble`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
DGEdat_PCcor <- cbind(DGEdat.cpm$samples$lib.size, DGEdat_PCcor)
rownames(DGEdat_PCcor) <- DGEdat.cpm$samples$sample_id

# Make a dataframe containing all plotting info
percentVar <- pca$sdev^2/sum(pca$sdev^2) # Save % variation by PC1 and PC2
d <- data.frame(PC1 = pca$x[, 1], 
                PC2 = pca$x[, 2],
                Group = DGEdat.cpm$samples[[group]],
                DGEdat.cpm$samples, name = colnames(DGEdat.cpm$counts)
              )

allgenes_PCA <- ggplot(data = d, aes_string(x = "PC1", y = "PC2", colour = "Group")) + # NOTE: can add 'shape' to be appropriate metadata column to add extra info
  geom_point(size = 3) + 
  xlab(paste0("PC1: ", round(percentVar[1] * 100), "% variance")) + 
  ylab(paste0("PC2: ", round(percentVar[2] * 100), "% variance")) + 
  coord_fixed() + 
  theme_bw() + # Set background color
  theme(panel.border = element_blank(), # Set border
        panel.grid.major = element_blank(), # Set major gridlines
        panel.grid.minor = element_blank(), # Set minor gridlines
        axis.line = element_line(colour = "black", size = 0.6), # Set axes color
        plot.background=element_blank(), # Set the plot background
        axis.title = element_text(size = 14), # Axis title size
        axis.text = element_blank()) # Axis text size and view plot
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
allgenes_PCA

6. Compile WGCNA Dataset

Transpose the filtered gene count matrix so that the gene IDs are rows and the sample IDs are columns.

datExpr <- as.data.frame(t(DGEdat.cpm$counts))

Check for genes and samples with too many missing values with goodSamplesGenes. There shouldn’t be any because we performed pre-filtering

gsg = goodSamplesGenes(datExpr, verbose = 3)
##  Flagging genes and samples with too many missing values...
##   ..step 1
gsg$allOK # Should return TRUE if not, the R chunk below will take care of flagged data
## [1] TRUE

Remove flagged samples if the allOK is FALSE

print(paste("Number genes before filtering: ", ncol(datExpr), sep='')) # Number genes before
## [1] "Number genes before filtering: 15467"
if (!gsg$allOK) # If the allOK is FALSE...
{
  # Optionally, print the gene and sample names that are flagged:
  if (sum(!gsg$goodGenes)>0) {
    printFlush(paste("Removing genes:", paste(names(datExpr)[!gsg$goodGenes], collapse = ", ")))
  }
  if (sum(!gsg$goodSamples)>0) {
    printFlush(paste("Removing samples:", paste(rownames(datExpr)[!gsg$goodSamples], collapse = ", ")))
  }
  # Remove the offending genes and samples from the data:
  datExpr = datExpr[gsg$goodSamples, gsg$goodGenes]
}
print(paste("Number genes after filtering:  ", ncol(datExpr), sep='')) # Number genes after
## [1] "Number genes after filtering:  15467"

7. Cluster the samples to look for obvious outliers

Look for outliers by examining the following dendrogram.

sampleTree = hclust(dist(datExpr), method = "average");
# Plot the sample tree: Open a graphic output window of size 12 by 8 inches
# The user should change the dimensions if the window is too large or too small.
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5, cex.axis = 1.5, cex.main = 2)

Exclude selected samples. We will exlcude Pacuta_HTHC_TP11_1416 and Pacuta_HTAC_TP9_1302 because they appear to be outliers on the dendrogram.

print(paste("Number samples before outlier removal: ", nrow(datExpr), sep='')) # Number samples before
## [1] "Number samples before outlier removal: 119"
row_names_df_to_remove<-c("Pacuta_HTHC_TP11_1416", "Pacuta_HTAC_TP9_1302") # Picked from above plot
datExpr <- (datExpr[!(row.names(datExpr) %in% row_names_df_to_remove),])
treatmentinfo <- (treatmentinfo[!(treatmentinfo$sample_id %in% row_names_df_to_remove),])

print(paste("Number samples after outlier removal:  ", nrow(datExpr), sep='')) # Number samples after
## [1] "Number samples after outlier removal:  117"

8. Network construction and consensus module detection

Choosing a soft-thresholding power: Analysis of network topology β

The soft thresholding power (β) is the number to which the co-expression similarity is raised to calculate adjacency. The function pickSoftThreshold performs a network topology analysis. The user chooses a set of candidate powers, however the default parameters are suitable values.

# Choose a set of soft-thresholding powers
powers <- c(seq(from = 1, to=10, by=0.5)) # Create a string of numbers from 1 through 10, increasing my 0.5 incraments

# Call the network topology analysis function
sft <-pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)
## pickSoftThreshold: will use block size 2892.
##  pickSoftThreshold: calculating connectivity for given powers...
##    ..working on genes 1 through 2892 of 15467
## Warning: executing %dopar% sequentially: no parallel backend registered
##    ..working on genes 2893 through 5784 of 15467
##    ..working on genes 5785 through 8676 of 15467
##    ..working on genes 8677 through 11568 of 15467
##    ..working on genes 11569 through 14460 of 15467
##    ..working on genes 14461 through 15467 of 15467
##    Power SFT.R.sq  slope truncated.R.sq mean.k. median.k. max.k.
## 1    1.0   0.0345  0.574          0.970 2820.00   2750.00   4700
## 2    1.5   0.0626 -0.534          0.950 1460.00   1390.00   3100
## 3    2.0   0.3460 -1.110          0.951  817.00    747.00   2170
## 4    2.5   0.6020 -1.460          0.964  486.00    423.00   1600
## 5    3.0   0.7330 -1.710          0.962  304.00    249.00   1210
## 6    3.5   0.8100 -1.830          0.969  198.00    152.00    944
## 7    4.0   0.8440 -1.900          0.966  134.00     95.30    750
## 8    4.5   0.8610 -1.940          0.966   92.70     61.00    605
## 9    5.0   0.8820 -1.940          0.971   66.10     40.00    495
## 10   5.5   0.8910 -1.930          0.974   48.10     26.60    409
## 11   6.0   0.9000 -1.930          0.978   35.80     18.10    346
## 12   6.5   0.9050 -1.920          0.982   27.10     12.50    294
## 13   7.0   0.9090 -1.900          0.985   20.80      8.71    253
## 14   7.5   0.9080 -1.890          0.984   16.20      6.22    218
## 15   8.0   0.9030 -1.880          0.983   12.80      4.48    190
## 16   8.5   0.9030 -1.860          0.985   10.20      3.25    166
## 17   9.0   0.8990 -1.860          0.984    8.28      2.39    146
## 18   9.5   0.8940 -1.840          0.983    6.75      1.77    129
## 19  10.0   0.8940 -1.820          0.984    5.56      1.32    114

Plot the results.

sizeGrWindow(9, 5)
par(mfrow = c(1,2));
cex1 = 0.9;
# Scale-free topology fit index as a function of the soft-thresholding power
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     xlab="Soft Threshold (power)",
     ylab="Scale Free Topology Model Fit,signed R^2",
     type="n",
     main = paste("Scale independence")
    )
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
    labels=powers,
    cex=cex1,
    col="red"
  )
# This line corresponds to using an R^2 cut-off
abline(h=0.8,col="red")
# Mean connectivity as a function of the soft-thresholding power
plot(sft$fitIndices[,1], sft$fitIndices[,5],
    xlab="Soft Threshold (power)",
    ylab="Mean Connectivity",
    type="n",
    main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5],
     labels=powers,
     cex=cex1,
     col="red"
    )

Select best soft power threshold using the above plot. I used a scale-free topology fit index R^2 of 0.831. This lowest recommended R^2 by Langfelder and Horvath is 0.8. I chose 0.83 because we want to use the smallest soft thresholding power that maximizes with model fit. It appears that our soft thresholding power is 3.5 because it is the lowest power before the R^2=0.8 threshold that maximizes with model fit.

Set softPower

softPower <- 3.5

9. Step-by-step network construction and module detection:

Co-expression adjacency and topological overlap matrix similarity

The next few steps may need to be executed on a supercomputer as our dataset is too large for most standard laptops to handle.

# Save data (on local machine)
#save(datExpr, file = "2a-WGCNA.RData")

# Load data (on server)
#adjTOM <- load(file="2a-WGCNA.RData")

May have to execute an the HPC server
Co-expression similarity and adjacency, using the soft thresholding power saved in softPower and translate the adjacency into topological overlap matrix to calculate the corresponding dissimilarity. I will use a signed network. https://peterlangfelder.com/2018/11/25/__trashed/

adjacency <- adjacency(datExpr, power=softPower,type="signed") # Calculate adjacency
TOM <- TOMsimilarity(adjacency, TOMType = "signed") # Translate adjacency into topological overlap matrix
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.
dissTOM <- 1-TOM # Calculate dissimilarity in TOM

Clustering using TOM

Form distance matrix.

geneTree <- flashClust(as.dist(dissTOM), method="average")

We will now plot a dendrogram of genes. Each leaf corresponds to a gene, branches grouping together densely are interconnected, highly co-expressed genes.

plot(geneTree, xlab="", sub="", main="Gene Clustering on TOM-based dissimilarity", labels= FALSE, hang=0.04)

Module identification using dynamicTreeCut

Module identification is essentially cutting the branches off the tree in the dendrogram above. We like large modules, so we set the minimum module size (minModuleSize at the top of script) relatively high. 30 is a good starting point. Module 0 is reserved for unassigned genes. The are other modules will be listed largest to smallest.

dynamicMods <- cutreeDynamic(dendro = geneTree,
                             distM = dissTOM,
                             deepSplit = 2,
                             pamRespectsDendro = FALSE,
                             minClusterSize = minModuleSize
                            )
##  ..cutHeight not given, setting it to 0.891  ===>  99% of the (truncated) height range in dendro.
##  ..done.
table(dynamicMods) # List modules and respective sizes
## dynamicMods
##    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16 
## 4315 3886 1409  869  673  529  449  365  338  333  311  237  214  203  175  170 
##   17   18   19   20   21   22   23   24   25 
##  146  142  139  136  134   98   74   66   56

Save results and reload on local machine if needed.

# Save data (on server)
#save(dynamicMods, geneTree, file = "2a-WGCNA.dyMod_geneTree.RData")

# Load data (on local machine)
#dyMod_geneTree <- load(file = "2a-WGCNA.dyMod_geneTree.RData")

Plot the module assignment under the gene dendrogram.

dynamicColors <- labels2colors(dynamicMods) # Convert numeric labels into colors
table(dynamicColors)
## dynamicColors
##         black          blue         brown          cyan     darkgreen 
##           449          3886          1409           203            98 
##      darkgrey       darkred darkturquoise         green   greenyellow 
##            66           134            74           673           311 
##        grey60     lightcyan    lightgreen   lightyellow       magenta 
##           146           170           142           139           338 
##  midnightblue        orange          pink        purple           red 
##           175            56           365           333           529 
##     royalblue        salmon           tan     turquoise        yellow 
##           136           214           237          4315           869
plotDendroAndColors(geneTree, dynamicColors, "Dynamic Tree Cut", dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05, main = "Gene dendrogram and module colors")

Merge modules whose expression profiles are very similar or choose not to merge

Plot module similarity based on eigengene value.

# Calculate eigengenes
MEList = moduleEigengenes(datExpr, colors = dynamicColors, softPower = 4.5)
MEs = MEList$eigengenes

# Calculate dissimilarity of module eigengenes
MEDiss = 1-cor(MEs)

# Cluster again and plot the results
METree = flashClust(as.dist(MEDiss), method = "average")

plot(METree, main = "Clustering of module eigengenes", xlab = "", sub = "")

ncol(MEs)
## [1] 25

10. Relating modules to developmental stage

Quantifying module–trait associations

Prepare trait data. Data has to be numeric, so I will substitute time_points and type for numeric values.

treatmentinfo$num <- c("1")
allTraits <- as.data.frame(pivot_wider(treatmentinfo, names_from = !!rlang::sym(group), values_from = num, id_cols = sample_id))
allTraits[is.na(allTraits)] <- c("0")
datTraits <- allTraits %>% column_to_rownames("sample_id")

Define numbers of genes and samples.

nGenes = ncol(datExpr)
nSamples = nrow(datExpr)

Correlations of traits with eigengenes.

moduleTraitCor = cor(MEs, datTraits, use = "p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples);
Colors=sub("ME","",names(MEs))

moduleTraitTree = hclust(dist(t(moduleTraitCor)), method = "average");
plot(moduleTraitTree, main = paste("'",group,"'"," clustering based on module-trait correlation", sep=''), sub="", xlab="", cex.lab = 1.5, cex.axis = 1.5, cex.main = 2)

Correlations of genes with eigengenes.

moduleGeneCor <- cor(MEs,datExpr)
moduleGenePvalue <- corPvalueStudent(moduleGeneCor, nSamples)

Plot module-trait associations

Represent module trait correlations as a heatmap.

textMatrix <- paste(signif(moduleTraitCor, 2), "\n(",signif(moduleTraitPvalue, 1), ")", sep = "")
dim(textMatrix) <- dim(moduleTraitCor)
head(textMatrix)
##      [,1]            [,2]           [,3]             [,4]            
## [1,] "-0.048\n(0.6)" "0.14\n(0.1)"  "0.047\n(0.6)"   "0.072\n(0.4)"  
## [2,] "-0.13\n(0.2)"  "0.066\n(0.5)" "0.15\n(0.1)"    "0.13\n(0.2)"   
## [3,] "-0.0056\n(1)"  "0.0022\n(1)"  "-0.15\n(0.1)"   "-0.25\n(0.006)"
## [4,] "-0.069\n(0.5)" "0.18\n(0.06)" "0.096\n(0.3)"   "0.018\n(0.8)"  
## [5,] "0.032\n(0.7)"  "0.11\n(0.2)"  "-0.66\n(4e-16)" "-0.14\n(0.1)"  
## [6,] "0.22\n(0.02)"  "0.16\n(0.09)" "-0.3\n(9e-04)"  "-0.43\n(1e-06)"
##      [,5]            [,6]            [,7]            [,8]            
## [1,] "-0.21\n(0.02)" "-0.08\n(0.4)"  "0.091\n(0.3)"  "-0.12\n(0.2)"  
## [2,] "-0.041\n(0.7)" "-0.092\n(0.3)" "0.033\n(0.7)"  "-0.33\n(3e-04)"
## [3,] "0.13\n(0.2)"   "0.55\n(2e-10)" "-0.045\n(0.6)" "-0.32\n(4e-04)"
## [4,] "-0.23\n(0.01)" "-0.039\n(0.7)" "0.086\n(0.4)"  "-0.15\n(0.1)"  
## [5,] "0.075\n(0.4)"  "0.69\n(7e-18)" "-0.039\n(0.7)" "-0.013\n(0.9)" 
## [6,] "0.12\n(0.2)"   "0.56\n(4e-11)" "-0.065\n(0.5)" "-0.081\n(0.4)" 
##      [,9]           
## [1,] "-0.047\n(0.6)"
## [2,] "0.087\n(0.4)" 
## [3,] "0.14\n(0.1)"  
## [4,] "-0.081\n(0.4)"
## [5,] "0.081\n(0.4)" 
## [6,] "0.043\n(0.6)"
labeledHeatmap(Matrix = moduleTraitCor,
               xLabels = names(datTraits),
               yLabels = names(MEs),
               ySymbols = names(MEs),
               cex.lab.y= 0.55,
               cex.lab.x= 0.55,
               colors = blueWhiteRed(50),
               textMatrix = textMatrix,
               setStdMargins = TRUE,
               cex.text = 0.25,
               textAdj = 0,
               zlim = c(-1,1),
               main = paste("Module-trait relationships")
              )

Represent module trait correlations as a complexHeatmap

# Create list of p-values for eigengene correlation with specific life stages
heatmappval <- signif(moduleTraitPvalue, 1)

# Make list of heatmap row colors
htmap.colors <- names(MEs)
htmap.colors <- gsub("ME", "", htmap.colors)

ht=Heatmap(moduleTraitCor,
           name = "Eigengene",
           column_title = "Module-Trait Eigengene Correlation",
           col = blueWhiteRed(50),
           row_names_side = "left", row_dend_side = "left",
           width = unit(5, "in"), height = unit(8.5, "in"),
           column_order = 1:ncol(moduleTraitCor),
           column_dend_reorder = TRUE,
           cluster_columns = hclust(dist(t(moduleTraitCor)), method = "average"),
           column_dend_height = unit(0.5, "in"),
           cluster_rows = METree, row_gap = unit(2.5, "mm"), border = TRUE,
           cell_fun = function(j, i, x, y, w, h, col) {
           if(heatmappval[i, j] <= 0.05) {
               grid.text(sprintf("%s", heatmappval[i, j]), x, y, gp = gpar(fontsize = 8, fontface = "bold"))
           } else {
               grid.text(sprintf("%s", heatmappval[i, j]), x, y, gp = gpar(fontsize = 8, fontface = "plain"))
           }},
           column_names_gp =  gpar(fontsize = 10),
           row_names_gp = gpar(fontsize = 10, alpha = 0.75, border = TRUE, fill = htmap.colors)
          )
draw(ht)

11. Gene relationship to trait and important modules: Gene Significance and Module Membership

We quantify associations of individual genes with life stage by defining Gene Significance GS as the absolute value of the correlation between the gene and the trait in group. For each module, we also define a quantitative measure of module membership MM as the correlation of the module eigengene and the gene expression profile.

Define variable weight containing the weight column of datTrait.

trait <- as.data.frame(as.numeric(as.factor(treatmentinfo[[group]]))) 
names(trait) = "group"
dim(trait)
## [1] 117   1

Colors of the modules.

modNames <- substring(names(MEs), 3)

geneModuleMembership <- as.data.frame(cor(datExpr, MEs, use = "p"))
MMPvalue <- as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples))

names(geneModuleMembership) <- paste("MM", modNames, sep="")
names(MMPvalue) <- paste("p.MM", modNames, sep="")

geneTraitSignificance <- as.data.frame(cor(datExpr, trait, use = "p"))
GSPvalue <- as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples))

names(geneTraitSignificance) <- paste("GS.", names(trait), sep="")
names(GSPvalue) <- paste("p.GS.", names(trait), sep="")

12. Summary output of network analysis results

Create the starting data frame.

geneInfo0 = data.frame(
  gene_id = names(datExpr),
  moduleColor = dynamicColors,
  geneTraitSignificance,
  GSPvalue
)

Order modules by their significance for trait.

modOrder = order(-abs(cor(MEs, trait, use = "p")))

Add module membership information in the chosen order.

for (mod in 1:ncol(geneModuleMembership))
{
  oldNames = names(geneInfo0)
  geneInfo0 = data.frame(geneInfo0, 
                         geneModuleMembership[, modOrder[mod]],
                         MMPvalue[, modOrder[mod]]
                         )
  names(geneInfo0) = c(oldNames, 
                       paste("MM.", modNames[modOrder[mod]], sep=""),
                       paste("p.MM.", modNames[modOrder[mod]], sep="")
                       )
}

Order the genes in the geneInfo variable first by module color, then by geneTraitSignificance.

geneOrder = order(geneInfo0$moduleColor, -abs(geneInfo0[[paste("p.GS.", names(trait), sep="")]]));
geneInfo = geneInfo0[geneOrder, ]
head(geneInfo)
##                                                                                 gene_id
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1     Pocillopora_acuta_KBHIv2___RNAseq.g880.t1
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1   Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t       Pocillopora_acuta_KBHIv2___RNAseq.9444_t
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1
##                                             moduleColor      GS.group
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1       black -2.068701e-05
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1         black -6.102154e-04
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1        black  1.806745e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1       black  2.002178e-03
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t          black -2.866394e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1       black  3.701016e-03
##                                             p.GS.group     MM.green p.MM.green
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1  0.9998234  0.143501397 0.12269986
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1    0.9947902 -0.011848636 0.89910582
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1   0.9845754  0.155248958 0.09464378
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1  0.9829071 -0.080992897 0.38534321
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t     0.9755312 -0.010973204 0.90652525
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1  0.9684097 -0.004831206 0.95877036
##                                               MM.yellow  p.MM.yellow
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.34892453 0.0001154571
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1   -0.30995655 0.0006717571
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1  -0.04759341 0.6103575539
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.08257903 0.3760748770
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t    -0.03230084 0.7295486759
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.03111254 0.7391335203
##                                              MM.magenta p.MM.magenta
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.06055089  0.516647736
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1   -0.13133985  0.158089258
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1   0.28681871  0.001717039
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.07031045  0.451274908
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t    -0.19967788  0.030893385
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.02746369  0.768810515
##                                             MM.lightcyan p.MM.lightcyan
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1  -0.09183281    0.324751101
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1    -0.04236013    0.650203857
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1   -0.25115627    0.006306037
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1  -0.09893561    0.288567834
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t      0.11024440    0.236697001
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1   0.08171378    0.381113873
##                                             MM.darkgrey p.MM.darkgrey
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.39969467  8.026113e-06
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1   -0.25671100  5.207297e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1  -0.01807584  8.466170e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.09963858  2.851386e-01
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t    -0.02499447  7.890876e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.09713966  2.974529e-01
##                                                MM.brown   p.MM.brown
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.48062039 4.150939e-08
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1   -0.50966309 4.389294e-09
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1  -0.11151680 2.312966e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.17399406 6.063274e-02
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t    -0.02069013 8.247663e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.12194364 1.902791e-01
##                                             MM.darkgreen p.MM.darkgreen
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1  -0.20760835     0.02470093
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1    -0.09794059     0.29346839
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1    0.07561684     0.41776749
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1   0.04357657     0.64084788
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t     -0.01772474     0.84956079
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1  -0.01279809     0.89106891
##                                             MM.royalblue p.MM.royalblue
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1  -0.07197915    0.440579294
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1    -0.13942809    0.133799310
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1   -0.29455654    0.001265272
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1  -0.09559146    0.305255539
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t      0.15729986    0.090312512
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1   0.05430192    0.560911559
##                                             MM.darkturquoise p.MM.darkturquoise
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1      -0.33623611       2.102125e-04
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1        -0.37199972       3.630637e-05
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1        0.01426237       8.786959e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1      -0.30247252       9.177075e-04
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t         -0.24991345       6.578395e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1      -0.13125934       1.583466e-01
##                                              MM.darkred p.MM.darkred
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.38741657 1.593677e-05
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1   -0.30741407 7.475506e-04
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1   0.02616069 7.794922e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.06924199 4.581987e-01
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t    -0.21558140 1.957728e-02
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.06917848 4.586121e-01
##                                             MM.midnightblue p.MM.midnightblue
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1     -0.57279046      1.498926e-11
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1       -0.48048974      4.191157e-08
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1      -0.04847061      6.037845e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1     -0.36559497      5.050166e-05
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t        -0.29965632      1.029850e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1     -0.23009610      1.256980e-02
##                                               MM.grey60  p.MM.grey60
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.49092055 1.916181e-08
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1   -0.57232177 1.570571e-11
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1  -0.05451036 5.594064e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.27574594 2.619439e-03
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t    -0.12113522 1.932563e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.17041330 6.621404e-02
##                                             MM.lightgreen p.MM.lightgreen
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1   -0.30367702    0.0008732472
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1     -0.32251816    0.0003905815
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1     0.18435207    0.0466183089
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1   -0.12317429    0.1858110291
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t      -0.26891663    0.0033705434
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1   -0.09830956    0.2916447722
##                                             MM.purple  p.MM.purple
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 0.6406916 7.281856e-15
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1   0.3125822 6.009407e-04
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1  0.2499209 6.576734e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 0.2765440 2.542343e-03
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t    0.1864207 4.416974e-02
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 0.3168960 4.993218e-04
##                                             MM.greenyellow p.MM.greenyellow
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1     -0.6903599     7.255620e-18
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1       -0.6298564     2.791317e-14
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1      -0.1537151     9.798913e-02
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1     -0.2641491     4.004189e-03
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t        -0.2392598     9.375404e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1     -0.2517965     6.169689e-03
##                                             MM.lightyellow p.MM.lightyellow
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1      0.6421866     6.024165e-15
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1        0.3666212     4.792289e-05
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1       0.2358117     1.048198e-02
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1      0.2110167     2.238515e-02
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t         0.1574486     9.000466e-02
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1      0.3705038     3.923999e-05
##                                                 MM.pink  p.MM.pink     MM.tan
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.20804558 0.02439285 -0.5528863
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1   -0.18357048 0.04757232 -0.4514506
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1   0.20585518 0.02596985 -0.1094448
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1  0.05946854 0.52418457 -0.1786946
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t    -0.17070584 0.06574293 -0.1504302
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.02806364 0.76390665 -0.1746129
##                                                 p.MM.tan    MM.red     p.MM.red
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 1.022145e-10 0.6079753 3.607703e-13
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1   3.244156e-07 0.4803911 4.221761e-08
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1  2.401354e-01 0.1739004 6.077368e-02
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 5.389617e-02 0.4355578 9.199729e-07
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t    1.054672e-01 0.1539702 9.742643e-02
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 5.970825e-02 0.2706109 3.168059e-03
##                                             MM.turquoise p.MM.turquoise
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1    0.3097140   0.0006786718
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1      0.3345102   0.0002276165
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1     0.2110974   0.0223326587
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1    0.2348850   0.0107982560
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t       0.1835526   0.0475943687
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1    0.1292410   0.1649004158
##                                               MM.orange p.MM.orange     MM.blue
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.18192643 0.049631983  0.07822197
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1   -0.25421504 0.005677917 -0.17397229
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1   0.16746756 0.071112581  0.12847846
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1  0.18236747 0.049072337  0.06880252
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t    -0.02848361 0.760479455  0.07526046
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.14915004 0.108499564  0.09219025
##                                              p.MM.blue  MM.black   p.MM.black
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 0.40186157 0.7200100 5.803009e-20
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1   0.06066546 0.5701121 1.955352e-11
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1  0.16742817 0.2922912 1.384786e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 0.46106352 0.5640874 3.524952e-11
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t    0.41997154 0.2995879 1.032722e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 0.32286342 0.3053941 8.132671e-04
##                                               MM.cyan    p.MM.cyan   MM.salmon
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 0.4415957 6.230095e-07 -0.60650930
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1   0.3530400 9.453814e-05 -0.68246587
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1  0.2194371 1.744741e-02 -0.08393264
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 0.5526485 1.045064e-10 -0.36854466
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t    0.2462552 7.441642e-03 -0.16963764
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 0.1505055 1.052907e-01 -0.11927007
##                                              p.MM.salmon
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 4.252695e-13
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1   2.383775e-17
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1  3.682737e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 4.341841e-05
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t    6.747649e-02
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 2.002537e-01

Show and save geneInfo as a TSV.

dim(geneInfo)
## [1] 15467    54
head(geneInfo)
##                                                                                 gene_id
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1     Pocillopora_acuta_KBHIv2___RNAseq.g880.t1
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1   Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t       Pocillopora_acuta_KBHIv2___RNAseq.9444_t
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1
##                                             moduleColor      GS.group
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1       black -2.068701e-05
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1         black -6.102154e-04
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1        black  1.806745e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1       black  2.002178e-03
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t          black -2.866394e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1       black  3.701016e-03
##                                             p.GS.group     MM.green p.MM.green
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1  0.9998234  0.143501397 0.12269986
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1    0.9947902 -0.011848636 0.89910582
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1   0.9845754  0.155248958 0.09464378
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1  0.9829071 -0.080992897 0.38534321
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t     0.9755312 -0.010973204 0.90652525
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1  0.9684097 -0.004831206 0.95877036
##                                               MM.yellow  p.MM.yellow
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.34892453 0.0001154571
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1   -0.30995655 0.0006717571
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1  -0.04759341 0.6103575539
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.08257903 0.3760748770
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t    -0.03230084 0.7295486759
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.03111254 0.7391335203
##                                              MM.magenta p.MM.magenta
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.06055089  0.516647736
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1   -0.13133985  0.158089258
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1   0.28681871  0.001717039
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.07031045  0.451274908
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t    -0.19967788  0.030893385
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.02746369  0.768810515
##                                             MM.lightcyan p.MM.lightcyan
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1  -0.09183281    0.324751101
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1    -0.04236013    0.650203857
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1   -0.25115627    0.006306037
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1  -0.09893561    0.288567834
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t      0.11024440    0.236697001
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1   0.08171378    0.381113873
##                                             MM.darkgrey p.MM.darkgrey
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.39969467  8.026113e-06
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1   -0.25671100  5.207297e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1  -0.01807584  8.466170e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.09963858  2.851386e-01
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t    -0.02499447  7.890876e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.09713966  2.974529e-01
##                                                MM.brown   p.MM.brown
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.48062039 4.150939e-08
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1   -0.50966309 4.389294e-09
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1  -0.11151680 2.312966e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.17399406 6.063274e-02
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t    -0.02069013 8.247663e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.12194364 1.902791e-01
##                                             MM.darkgreen p.MM.darkgreen
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1  -0.20760835     0.02470093
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1    -0.09794059     0.29346839
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1    0.07561684     0.41776749
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1   0.04357657     0.64084788
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t     -0.01772474     0.84956079
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1  -0.01279809     0.89106891
##                                             MM.royalblue p.MM.royalblue
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1  -0.07197915    0.440579294
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1    -0.13942809    0.133799310
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1   -0.29455654    0.001265272
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1  -0.09559146    0.305255539
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t      0.15729986    0.090312512
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1   0.05430192    0.560911559
##                                             MM.darkturquoise p.MM.darkturquoise
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1      -0.33623611       2.102125e-04
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1        -0.37199972       3.630637e-05
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1        0.01426237       8.786959e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1      -0.30247252       9.177075e-04
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t         -0.24991345       6.578395e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1      -0.13125934       1.583466e-01
##                                              MM.darkred p.MM.darkred
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.38741657 1.593677e-05
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1   -0.30741407 7.475506e-04
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1   0.02616069 7.794922e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.06924199 4.581987e-01
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t    -0.21558140 1.957728e-02
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.06917848 4.586121e-01
##                                             MM.midnightblue p.MM.midnightblue
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1     -0.57279046      1.498926e-11
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1       -0.48048974      4.191157e-08
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1      -0.04847061      6.037845e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1     -0.36559497      5.050166e-05
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t        -0.29965632      1.029850e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1     -0.23009610      1.256980e-02
##                                               MM.grey60  p.MM.grey60
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.49092055 1.916181e-08
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1   -0.57232177 1.570571e-11
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1  -0.05451036 5.594064e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 -0.27574594 2.619439e-03
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t    -0.12113522 1.932563e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.17041330 6.621404e-02
##                                             MM.lightgreen p.MM.lightgreen
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1   -0.30367702    0.0008732472
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1     -0.32251816    0.0003905815
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1     0.18435207    0.0466183089
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1   -0.12317429    0.1858110291
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t      -0.26891663    0.0033705434
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1   -0.09830956    0.2916447722
##                                             MM.purple  p.MM.purple
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 0.6406916 7.281856e-15
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1   0.3125822 6.009407e-04
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1  0.2499209 6.576734e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 0.2765440 2.542343e-03
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t    0.1864207 4.416974e-02
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 0.3168960 4.993218e-04
##                                             MM.greenyellow p.MM.greenyellow
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1     -0.6903599     7.255620e-18
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1       -0.6298564     2.791317e-14
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1      -0.1537151     9.798913e-02
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1     -0.2641491     4.004189e-03
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t        -0.2392598     9.375404e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1     -0.2517965     6.169689e-03
##                                             MM.lightyellow p.MM.lightyellow
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1      0.6421866     6.024165e-15
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1        0.3666212     4.792289e-05
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1       0.2358117     1.048198e-02
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1      0.2110167     2.238515e-02
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t         0.1574486     9.000466e-02
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1      0.3705038     3.923999e-05
##                                                 MM.pink  p.MM.pink     MM.tan
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.20804558 0.02439285 -0.5528863
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1   -0.18357048 0.04757232 -0.4514506
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1   0.20585518 0.02596985 -0.1094448
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1  0.05946854 0.52418457 -0.1786946
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t    -0.17070584 0.06574293 -0.1504302
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.02806364 0.76390665 -0.1746129
##                                                 p.MM.tan    MM.red     p.MM.red
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 1.022145e-10 0.6079753 3.607703e-13
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1   3.244156e-07 0.4803911 4.221761e-08
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1  2.401354e-01 0.1739004 6.077368e-02
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 5.389617e-02 0.4355578 9.199729e-07
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t    1.054672e-01 0.1539702 9.742643e-02
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 5.970825e-02 0.2706109 3.168059e-03
##                                             MM.turquoise p.MM.turquoise
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1    0.3097140   0.0006786718
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1      0.3345102   0.0002276165
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1     0.2110974   0.0223326587
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1    0.2348850   0.0107982560
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t       0.1835526   0.0475943687
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1    0.1292410   0.1649004158
##                                               MM.orange p.MM.orange     MM.blue
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 -0.18192643 0.049631983  0.07822197
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1   -0.25421504 0.005677917 -0.17397229
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1   0.16746756 0.071112581  0.12847846
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1  0.18236747 0.049072337  0.06880252
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t    -0.02848361 0.760479455  0.07526046
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 -0.14915004 0.108499564  0.09219025
##                                              p.MM.blue  MM.black   p.MM.black
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 0.40186157 0.7200100 5.803009e-20
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1   0.06066546 0.5701121 1.955352e-11
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1  0.16742817 0.2922912 1.384786e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 0.46106352 0.5640874 3.524952e-11
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t    0.41997154 0.2995879 1.032722e-03
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 0.32286342 0.3053941 8.132671e-04
##                                               MM.cyan    p.MM.cyan   MM.salmon
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 0.4415957 6.230095e-07 -0.60650930
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1   0.3530400 9.453814e-05 -0.68246587
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1  0.2194371 1.744741e-02 -0.08393264
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 0.5526485 1.045064e-10 -0.36854466
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t    0.2462552 7.441642e-03 -0.16963764
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 0.1505055 1.052907e-01 -0.11927007
##                                              p.MM.salmon
## Pocillopora_acuta_KBHIv2___RNAseq.g13463.t1 4.252695e-13
## Pocillopora_acuta_KBHIv2___RNAseq.g880.t1   2.383775e-17
## Pocillopora_acuta_KBHIv2___RNAseq.g9596.t1  3.682737e-01
## Pocillopora_acuta_KBHIv2___RNAseq.g18377.t1 4.341841e-05
## Pocillopora_acuta_KBHIv2___RNAseq.9444_t    6.747649e-02
## Pocillopora_acuta_KBHIv2___RNAseq.g11595.t1 2.002537e-01
write.table(geneInfo, file=WGCNA_results.file, sep='\t', quote=F, row.names=F)
write.table(treatmentinfo, file=WGCNA_samples.file, sep='\t', quote=F, row.names=F)